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Abstract 

We present a methodology for determining reservoir groundwater age and tran- 
sit time probability distributions in a deterministic manner, considering advective- 
dispersive transport in steady velocity fields. In a first step, we propose to model 
the statistical distribution of groundwater age at aquifer scale by means of the 
classical advection-dispersion equation for a conservative and non-reactive tracer, 
associated to proper boundary conditions. The evaluated function corresponds to 
the density of probability of the random variable age, age being defined as the time 
elapsed since the water particles entered the aquifer. An adjoint backward model 
is introduced to characterize the life expectancy distribution, life expectancy being 
the time remaining before leaving the aquifer. By convolution of these two distribu- 
tions, groundwater transit time distributions, from inlet to outlet, are fully defined 
for the entire aquifer domain. In a second step, an accurate and efficient method 
is introduced to simulate the transit time distribution at discharge zones. By ap- 
plying the reservoir theory to advective-dispersive aquifer systems, we demonstrate 
that the discharge zone transit time distribution can be evaluated if the internal 
age probability distribution is known. The reservoir theory also applies to internal 
life expectancy probabilities yielding the recharge zone life expectancy distribution. 
Internal groundwater volumes are finally identified with respect to age and transit 
time. One- and two-dimensional theoretical examples are presented to illustrate the 
proposed mathematical models, and make inferences on the effect of aquifer struc- 
ture and macro-dispersion on the distributions of age, life expectancy and transit 
time. 
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1 Introduction. 



The knowledge of groundwater age distributions is of prime interest in many 
environmental issues since they depend on the intrinsic characteristics of the 
overall transport properties of an aquifer and its sub-systems. For instance, an 
important fraction of young water within a water sample can often be taken 
as the signature of a reservoir with good turnover property On the opposite, 
a considerable component of old water may reflect a poorly recharged aquifer, 
and/or significant internal mixing processes. The impact of a contamination 
hazard on groundwater quality can be investigated using groundwater age, 
since the age distribution provides direct information on the time required for 
a water particle, or a conservative tracer, to reach any critical zone that is 
to be protected. The fate of a solute being transported in groundwater par- 
tially depends on the time spent by the water molecules during their migration 
throughout the aquifer system. Reactive transport of a specific substance is 
also linked to groundwater age; the time spent flowing through any kind of 
mineral heterogeneity being a conditional factor for the development of any 
potential reactions. With the age information, inferences can be made on the 
aquifer physical characteristics, as well as on the chemical transformations 
that contaminants may undergo. The age can also be of importance for esti- 
mating historical aspects related e.g. to the agricultural practices, or the land 
use of a particular region, which are expected to have lead to groundwater 
contamination. As it is the case throughout this work, groundwater age can 
be considered as a fully conservative tracer. Consequently, the worst scenario 
with regard to a contamination case is thus chosen by solving the age transport 
problem. 

Generally, a water sample shows a mixture of ages, which may range between 
several orders of magnitude, as a consequence of the reservoir geometrical 
and hydro-dispersive characteristics (spatial repartition of the hydraulic and 
transport parameters). Groundwater age must, therefore, be regarded as a 
statistical, or probabilistic distribution, rather than considering it as a single 
absolute or average value. The dating methods commonly provide an average 
value over the water sample for the age of groundwater after recharge, which 
in theory does not represent hard data. In fact, the mean of an unknown dis- 
tribution, here the distribution of ages, is not necessarily a reliable value for 
the most likely of this distribution. By far the most frequently used dating 
methods are based on the measurements of natural tracers, such as the iso- 
topes of radon, carbon or oxygen, and on the measurements of man-induced 
atmospheric concentrations of elements such as tritium (^H), helium (^He), 
chloride (^^Cl), krypton (^^Kr), or chlorofluorocarbons (CFG), which have in- 
creased steadily between the 1940s and the early 1990s. The most efficient 
methods for dating recent waters are the ones based on ^H/^He, ^^Kr, and 
CFG- measurements [63,55,25], which are known to provide age dates over a 
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period of 40 years with an accuracy of 20% or less [11]. CFC-based ages pro- 
vide a good resolution for groundwater with relatively young components [8,6], 
but they do not account for the time spent by water molecules within the un- 
saturated zone, such that in the case of deep water table aquifers the travel 
time duration within the unsaturated zone must, therefore, be estimated. Ab- 
solute dating techniques involve decay of radionuclides in groundwater, for 
which ^^C or '^^Cl are classical elements that are used for dating old ground- 
water, e.g. in deep and large sedimentary basins. Groundwater age is very 
often used to make inferences on aquifer parameters, groundwater recharge, 
flow paths and flow velocities. However, with many of these dating methods, 
the interpretation of ages is achieved by means of simple conceptual models 
and fitting the data. This involves significant simplifications of the fiow and 
transport processes, which may lead to erroneous interpretations about e.g. 
the past release of contaminants in aquifers. As will be discussed in this paper, 
mixing and dispersion are major factors, which can lead to unrepresentative 
mean age measurements. Dating methods remain however very useful for cali- 
brating numerical models [56,58], which attempt to simulate the flow patterns, 
the flow rates, and the distribution of ages in groundwater. 

To quantify the distribution of ages in aquifers, several types of mathematical 
models have been developed during the past decades of research in this field. 
Nevertheless, in many cases age distributions are more or less arbitrarily cho- 
sen and not deterministically calculated. Analytical lumped-parameter type 
models have been however extensively used in the simulation of environmental 
tracer data [45,47,74,9,57,1], such as isotopic data, which are commonly inter- 
preted with advective transit time models, although isotope transport does not 
necessarily undergo advective processes only. Specified transit time distribu- 
tions describing piston-fiow, exponential mixing, combined piston-exponential 
mixing, or dispersive mixing, are chosen to solve the inverse problem by fitting 
the model on measured tracer output data. This procedure calls for significant 
simplifications, which can often not be justified, such as for instance neglecting 
the reservoir structure, as well as the spatial variability of infiltration rates [10] 
and aquifer fiow and transport parameters. Amin and Campana [1] proposed 
to model the groundwater age mixing process by means of a three-parameter 
gamma function which accounts for various states of mixing ranging between 
no mixing (piston fiow model) and perfect mixing (exponential model). Ro- 
bust verifications of the applicability of lumped parameter models can hardly 
be found [36,44,28]. Transit time distributions are often obtained with nu- 
merical solutions by making the assumption of pure advective motion of the 
groundwater particles, the particle-tracking technique being the most popular 
one. Purely kinematic ages ignore the effect of dispersion and mixing on age 
transport [22,12], and often reveal to be ill posed in complex heterogeneous 
systems [71], for which the 3-D implementation is subject to severe technical 
problems. The importance of including age dilution processes such as disper- 
sion and matrix diffusion, when comparison is made between modeled and 
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measured ages, has been pointed out by many authors [66,45,46]. Moreover, 
particle-tracking does not allow calculation of transit time distributions since 
groundwater volumes are not associated to simulated ages. 

More elaborated quantitative approaches consider age as a mass that is trans- 
ported by groundwater through volume-averaged temporal moment equations 
[64,37,35,71], in which the product of groundwater age with its mass (age 
mass pA) is the conserved quantity. Harvey and Gorelick [37], and Varni and 
Carrera [71] derived a set of recursive temporal moment equations, which are 
sequentially solved in order to simulate the transit time distributions from the 
n calculated moments. According to Harvey and Gorelick [37], the first five 
moments that characterize the accumulated mass, the mean, the variance, the 
skewness and the kurtosis of a breakthrough curve, respectively, may provide 
sufficient information to summarize the entire distribution. Since many natu- 
ral systems reveal a multi-modality of the age distribution within the reservoir 
and at the discharge zones, and since the shape of this distribution is a priori 
unknown, an infinity of moments would therefore theoretically be required to 
construct the entire distribution. 

In the literature one can find many terms, which relate to a specific time 
spent by water molecules within the aquifer. Use will be made throughout 
this work of the notion of transit time as the total residence time of water 
molecules within the aquifer, i.e. the age of these molecules when they exit 
the aquifer. The notion of travel time is rather used to characterize the time 
spent to travel between two arbitrary locations in the aquifer. Travel time 
probabilities have been a subject of high interest in many studies characteriz- 
ing solute transport in sub-surface hydrology [16,18,19,38,39]. The travel time 
probability is commonly defined as the response function to an instantaneous 
unit flux impulse [21,39]. In their transfer function approach of contaminant 
transport through unsaturated soil units. Jury and Roth [39] model tracer 
breakthrough curves with one-dimensional travel time probability functions. 
Shapiro and Cvetkovic [60], Dagan [19], and Dagan and Nguyen [20] derived 
the forward travel time probability for a mass of solute by using the La- 
grangian concept of particle displacement in porous media. The derivation 
of forward and backward models for location and travel time probability has 
become a classical mathematical approach for contaminant transport charac- 
terization and prediction [69,43,50,51]. The spreading of a contaminant mass 
is analyzed by following the random motion of solute particles, and to do so, 
the advection-dispersion equation (ADE) is assimilated to the Fokker-Planck 
(or forward Kolmogorov) equation. The expected resident concentration of a 
conservative tracer is taken as the probability density function for the location 
of a particle, at any time after having entered the system. 

The aim of the present work is to provide a general theoretical framework to 
model complete groundwater age distributions at aquifer scale in a determin- 
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istic way. The concept of age variability is associated to the concept of random 
variables and their distributions by using classical elements of probability, al- 
lowing the introduction of mathematical definitions for age, life expectancy 
and transit time statistical distributions. Forward and backward ADEs for 
conservative tracers are used to simulate the above-mentioned distributions at 
aquifer scale. By manipulating the ADEs, the reservoir theory [27] is expressed 
in order to characterize recharge and discharge zones transit time distributions 
with refined accuracy. The proposed models are illustrated and discussed by 
means of analytical and numerical analysis of one- and two-dimensional the- 
oretical flow conflgurations. 



2 The 'ages' of groundwater as space-dependent random variables. 



In this section we present the models allowing the calculation of the statistical 
distribution of groundwater age, life expectancy and transit time in arbitrary 
aquifers. 



2. 1 Definitions. 



The characterization of groundwater molecules with respect to a travel time 
within an aquifer system is fully dependent on the spatial reference from which 
this time is "measured". Usually the groundwater age is defined as a relative 
quantity with respect to a starting location where age is assumed to be zero. 
Use will be made throughout this work of three variations of terminology. For 
a given spatial position in the reservoir, the age (A) relates to the time elapsed 
since the water molecules entered the system at the recharge limits, where age 
is zero. For the same spatial position, the life expectancy (E) is defined as the 
time required for the water molecules to reach an outlet limit of the system. 
Life expectancy is therefore zero at an outlet. The transit time (T) finally 
refers to the total time required by the same water molecules to migrate from 
an inlet zone (T = E) to an outlet zone {T — A). In a, REV, the three vari- 
ables A, E and T are random variables, characterized by probability density 
functions (pdf) that can be regarded as the statistical occurrence of water 
molecules with respect to time, which could be observed in a groundwater 
sample if any analytical procedure would allow such measurements. 
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2.2 Age probability. 



The typical heterogeneity of aquifer systems involves strongly varying flow ve- 
locity fields, with multi-scale coherence lengths. The spatial variability of ve- 
locity and transport parameters induces a spreading of the contaminant mass. 
The tensor of macro-dispersion in the classical ADE accounts for the uncer- 
tainty in the transport prediction induced by mixing. Various studies relate to 
how the ADE is limited by the impact of physical and chemical heterogeneities 
on solute transport, such that up-scaling is not always satisfying [17,67,34]. 
If such heterogeneities arc present at aquifer scale the transport parameters 
should be time-dependent, but this time dependency may be relaxed when the 
correlation scales of the transport parameter random fields are finite [34,49] . 
The ADE with time-independent parameters holds only when the solute par- 
ticles have enough time to be distributed by dispersion between the flow lines. 
Since we are interested in solving the age transport problem at aquifer scale, 
we make the assumption that the ADE with time-independent transport pa- 
rameters (the parameters have reached their asymptotic values) can model the 
evolutional transport of the groundwater age and life expectancy distributions 
under steady-flow conditions. The modelled process applies to conservative 
and non-reactive tracers. 

Let us consider an aquifer domain Q in the three-dimensional space, with 
hydraulic recharge boundary r_, discharge boundary r+, and impermeable 
boundary Tq, as illustrated in Fig. 1. The boundary r+ corresponds to the 
open boundary of the system, through which a free exit of age mass is ex- 
pected. With respect to the above-mentioned considerations on contaminant 
spreading, it is convenient to describe the groundwater sample age distribution 
as a random variable associated to a probabilistic model. The age probability 
distribution at a position x inQ can be evaluated by solving the ADE when a 
unit pulse of conservative tracer is uniformly applied on the recharge area r_ . 
The resulting breakthrough curve is the probabilistic age distribution [21,39]. 
Making use of this property, we propose to model groundwater age and life 
expectancy pdfs by forward and backward transient-state transport equations, 
under steady-state hydraulic conditions. The pre-solution of the velocity fleld 
is performed by the following steady-state groundwater flow equation: 



where q is the Darcy flux vector [LT~^], which is valid for ideal tracers, H 
is the hydraulic head [L], qi and go arc fluid source and sink terms [T^"*^], 
respectively, and K is the tensor of hydraulic conductivity [LT~^]. The age 
pdf is then obtained by solving the following forward boundary value problem: 



V • q = gi - go, q = -KVH 



in Q, 



(1) 



d(t)gA 
dt 



V • qgA + V • BVqa + qiS{t) - qoQA 



in Q 



(2a) 
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gA{x, 0) = gA{x, oo) = in Q 
3A{x,t)-n—{q-n)6{t) on r_ 
Ja{x, t) ■ n — on Fq 



(2b) 
(2c) 
(2d) 



where gA{x,t) is the transported age pdf [T~^], J^(a;, t) is the total age mass 
flux vector [LT~^], D is the tensor of macro-dispersion [L^T~^], x = (x,y,z) 
is the vector of cartesian coordinates [L], t is time [T], = (f){x) is porosity 
or mobile water content [— ], n is a normal outward unit vector, and 6{t) is 
the time-Dirac delta function [T~^], which ensures a pure flux impulse on r_. 
The source term qiS{t) is meant for simulating a potential internal production 
of water (3-D) or 2-D horizontal aquifer conflgurations with an areal recharge 
intensity qi. The sink term QoOa may result from any internal extraction of 
groundwater. The tensor of macro-dispersion D = 0D* = D(cc) in Eq. (2a) is 
deflned by Bear [3]: 

q (g) q 

D = (aj, - Q;t) II II + Q;r||q||I + 0-DmI (3) 



where and are the longitudinal and transversal coefficients of dispersiv- 
ity [L], respectively, is the coefficient of molecular diffusion [L^T~^], and 
I is the identity matrix. The total age mass flux vector 3A{x,t) is classically 
deflned by the sum of the convective and dispersive fluxes: 

Ja{x, t) = qgA{x, t) - TfVgAix, t) (4) 



2.3 Life expectancy probability. 



The hfe expectancy probabihty distribution satisfles the adjoint backward 
model of Eq. (2a): 

= V • qgE + V • BVgE - QWe in (5a) 

gE{x, 0) = gE{x, oo) = in Q (5b) 

JE{x,t) ■ n— —{q- n)S{t) on r+ (5c) 

- DV^E(a;, t)-n^O on Eq (5d) 

where gE{x,t) is the transported life expectancy pdf, and where the total life 
expectancy mass flux vector JE{x,t) is 

Je{x, t) = -qgE{x, t) - BVgsix, t) (6) 
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Eq. (5a) is the formal adjoint of Eq. (2a) [32,2], the so-called "backward-in- 
time" equation [69,73,72], also backward Kolmogorov equation [40]. Given the 
forward equation, the backward equation is technically obtained by reversing 
the sign of the flow fleld, and by adapting the boundary conditions [50,51,72]. 
On the impermeable boundary Fq, a third-kind condition (Cauchy) in the for- 
ward equation becomes a second-kind condition (Neumann) in the backward 
equation, and vice-versa. A second-kind condition in the forward model will 
also become a third-kind condition in the backward model [33, p. 146]. The 
advection term is known to be not self-adjoint (it should be written in the 
form q • Vqe in Eq. (5a)) unless flow is divergence free [72]. However, the 
backward equation can still handle non-divergence free flow fields by means 
of the important sink term —qiQE appearing in Eq. (5a). This sink term has 
been derived in Cornaton [13] from the vertical averaging process of the gen- 
eral 3-D backward ADE, and is consistent with the analysis of Neupauer and 
Wilson [51,52]. Recharge by internal sources (3-D or 2-D vertical) or by areal 
fluxes (fluid source for 2-D horizontal) is introduced by the first-order decay 
type term —qiQE, which is a consequence of the reversed fiow field. Internal 
sources produce a sink of life expectancy probability, while internal sinks (term 
QodA in Eq. (2a)) do not appear in the backward model since a fluid sink may 
not influence the hfe expectancy pdf. The boundary r_ corresponds to the 
open boundary of the system (since flow is reversed) through which a free exit 
of life expectancy mass is expected. 

The simulation of life expectancy with (5) is valid in the case of steady-state 
velocity fields only. Transient-state velocity fields would require another ap- 
propriate formulation of (5). For a steady-state hydraulic situation, if q ap- 
proaches zero in some regions of the reservoir, e.g. like in aquitards in which 
transport is diffusion-dominant, then (5) still holds because of the irreversible 
nature of dispersion. The amount of age mass diffused between an aquifer 
and an aquitard is proportional to concentration differences between the two 
media, and is the same in both the forward and backward problems. In the 
boundary value problems (2) and (5) the classical homogeneous Neumann 
boundary condition {—D'Vg ■ n — 0) at the outlet limit for the age problem 
(at the inlet limit for the life expectancy problem) is not used in order to 
allow a natural age/life expectancy gradient through the open boundaries. 
Instead, the normal projection of the dispersive ffux is evaluated implicitly 
at the boundaries. The evaluation procedure in the framework of the finite 
element method is described in Cornaton et al. [14]. This kind of boundary 
condition, which is referred to as Implicit Neumann condition, enables conti- 
nuity of the total mass flux at outlet. The Implicit Neumann condition is a 
generalized version of the Free Exit condition for parabolic equations proposed 
by Prind [31]. As discussed by Nauman and Buffham [48], Parker and van 
Genuchten [53], Kreft and Zuber [42], and Bear and Verruitj [4], total mass 
flux continuity at outlet permits upgradient solute movement by dispersion. 
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(a) (b) 

Fig. 1. Schematic illustration of a groundwater reservoir 0, with indicated no-flow 
(Fq), inlet (F_) and outlet (F_|_) boundaries: (a) Age problem with normal flow 
field; (b) Life expectancy problem with reversed flow fleld. The arrow heads on the 
symbolized flowlines (dashed lines) stand for the position and direction of water 
molecules at a given time after their release. The black dot stands for a small water 
sample, to illustrate the random variable transit time (T) as the sum of the two 
random variables age {A) and life expectancy (E). 

Eqs. (2) and (5) simulate the forward and backward transport resulting from a 
unit pulse input. The space-time evolution of the water molecules is described 
by the distributions qa^x, t) and gE{,x, t). Both differential equations deal with 
conditional probabilities that characterize the statistical occurrence of water 
molecules as a function of age and life expectancy. Location probability is 
related to resident concentration [41,18,39], and describes the position x of 
water molecules at a given time after their release in the system. On the other 
hand, travel time probability is related to flux concentration [41,18,39], and 
characterizes for a position x the amount of time spent within Q since the 
water molecules entered the system (in the right-hand side of Eq. (5c), age 
is at the flux concentration = S{t)). Resident concentration relates to the 
mass of solute per unit porous volume while flux concentration is defined as 
the solute mass flux per unit water flux. A flux concentration is the physical 
representation of the mean of the microscopic fluid concentrations weighted 
by the associated microscopic fluid velocities [53]. The multi-dimensional re- 
lation between flux and resident concentrations can be found in Sposito and 
Barry [65], and is formally the projection of the total mass flux on the flow 
velocity direction. Accordingly, the flux concentration form of the random 
variable age is: 

/ Ja-q Jl-q DV^A-q 

9 a = u ^ = 9A + J, ^ = Qa r, ^ (7) 

II q II II q II II q II 



with J = —J^Vg denoting the dispersive part of J. By analogy, the flux 
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(8) 



The age pdf at a position x inQ characterizes the probabihty per unit time for 
the time spent since recharge, and is a flux concentration, evaluated enforcing 
Eq. (7). The resident concentration of age characterizes a vohimc-avcragcd age 
density function, which describes the density of probabihty of finding water 
molecules at the position a; in at time t. 

2.4 Transit time probability. 

In Eqs. (2a) and (5a), the dependent variables are probability density functions 
of continuous time random variables. The behavior of these random variables 
can also be described by cumulative distribution functions (cdf). Let U be 
one of these two variables {U = A or U = E ), with u the associated values 
the variable U may take at a given position x of space. The cdf Gu{x, u) and 
the pdf gu{x,u) of the variable U are commonly defined as 



gu{x,u) = — ^ — , Gu{x, -00) = , Gu{x,oo) = 1 (10) 



where P denotes the probability event on U, or number of occurrences with 
U < u ratioed to the total number of occurrences, and r is a dummy vari- 
able for integration. The probability functions property (10) together with the 
boundary conditions in Eqs. (2) and (5) ensure that age and life expectancy 
qa and qe are directly related to probabilities, and since concentration can be 
modelled by the ADE, then probability too can be modelled by the ADE. For 
a given position x in ^2, the ages of groundwater molecules are described by 
the pdf Qai which measures the density of probability of having an age t. The 
same molecules are also characterized by the pdf ^^e, which measures the den- 
sity of probability of having a life expectancy t. Introducing now the random 
variable transit time T, with density of probability qt, the water molecules 
can be described by their intensity of probability of flowing throughout the 
system at a time t. The variable T is a random variable corresponding to the 
sum of the two random variables A and E (see Fig. 1). Hence the statistical 
distribution of T is the pdf of the sum of A and E, qt — Qa+e- This problem 




(9) 



with 
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can be solved if the joint pdf qa^e of A and E, which characterizes the joint 
behavior of A and is known [5] : 



The joint quantity gA,E{Xi a, ejdade relates to the probability that A lies in 
the small interval a to a + da, and that E lies in the small interval e to e + de. 
Assuming that A and E are stochastically independent variables (for the same 
spatial position, A depends on the initial point while E depends on the end 
point, and under a steady-state flow regime E may not be influenced by the 
memory of the past evolution), and since qa and qe are zero for negative values 
of their arguments, the joint pdf qa^e simplifies in QaOe [5]- The probability 
density function g-r in Eq. (11) can then be obtained using the convolution 
integral: 



from which the cdf of T can be calculated enforcing Eq. (9). The fact that 
qa and qe are zero for negative values of their arguments allows applying 
the convolution integral from to t. Since the pdfs qa and qe give the age 
and life expectancy probability of occurrence at each position x in Q, both 
the maximum age as well as the maximum life expectancy correspond to the 
maximum transit time. Consequently, the time variable t can equivalently refer 
to all specific values of age, life expectancy, and transit time. The convolution 
integral (12) states that the probabihty that the variable T lies in a small 
interval around t is proportional to the product of the probability that the 
variable A lies in the interval t to t + dt and a factor proportional to the 
probability that the variable E lies in a small interval around t — r, the value 
of E ensuring that A and E sum to T. This product is then summed over 
all possible values of time t (from the minimum age to the maximum age) to 
yield the transit time pdf at a position x in space. The derived distribution of 
T = A + E in Eq. (12) can also be viewed as a transfer function convolution 
process, the input distribution being the age pdf qa-, and the signal transferring 
function being the life expectancy qe- 

To our point of view, Eq. (12) is an important result of the present work. The 
field of Qt characterizes the evolution of groundwater molecules throughout the 
aquifer domain by specifying the amount of time from recharge to discharge. 
At a given position in the reservoir, the temporal evolution of the groundwater 
molecules can be characterized by the three pdfs qa, Qe and qt- Each function 
contains specific information on a time of residence, the nature of which is a 
function of the spatial references that are chosen for evaluation. For instance, 
qa is conditioned by the inlet limit r_, where the variable A is nil, while gs i^ 
conditioned by the outlet limit r+, where the variable E is nil. For the variable 




(11) 




(12) 
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Fig. 2. Age, life expectancy and transit time dimensionless pdfs in a 1-D domain 
for a Peclet number Pe = 20: (a) Flux pdfs; (b) Resident pdfs. Time is normalized 
by the average turnover time tq = L/v, x hy L, and Pe = Lv/D. The average age 
and the average life expectancy sum to the average transit time. 

T, the pdf qt is conditioned by the fact that T = A at outlet, and that T = E 
at inlet. 

If gA and qe are resident concentrations, so is gr- If Qa and qe are flux con- 
centrations, so is Qt- Applying a Laplace transform to Eqs. (7) and (8) and 
convoluting, the transit time resident and flux concentrations are found to be 
linked by the following relation: 

9^ = m = 9T-^, [gAi% - gE^A + ^^#^] ■ q (i3) 



where gu{x, s) denotes the s-transform state of the function gu{x,u), U = A, 
E or T. Eq. (13) shows that transit time flux concentration is dependent on the 
transit time resident concentration, but also on the tensor product between 
the age and life expectancy dispersive fluxes, and on the convolution product 
of the age and life expectancy flux and resident concentrations. For a zero 
dispersion case, = g for the random variables A, E and T. Consider the 
semi-infinite 1-D domain of characteristic length L (the outlet is supposed at 
the position x = L) and uniform velocity v along the illustrated in 

Fig. 2. The age flux pdf at the position x is obtained as the solution of Eq. (2a) 
using the boundary condition qa = gA = ^{t) at x = 0, and the age resident 
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pdf is obtained as the solution of Eq. (2a) using the boundary condition (2c) 
at a; = 0. These solutions are given in dimensionless form in Appendix A. 
They show the trivial fact that for a one- dimensional flow configuration, the 
transit time flux pdf is unique and independent on the spatial coordinate, as 
illustrated in Fig. 2a. For example, the age flux pdf at X = 1/4 in Fig. 2a is 
equal to the life expectancy flux pdf at X = 3/4, and the convolution of both 
distributions gives the transit time flux pdf, which equals both the age flux 
pdf at outlet (X = 1) and the life expectancy flux pdf at inlet (X = 0). The 
average age x/v and the average life expectancy {L — x)/v sum to the average 
transit time L/v, which is independent on x in a rectilinear flow line. In a 
similar way, the age resident pdf at the position X equals the life expectancy 
resident pdf at the position 1 — X, as shown in Fig. 2b. For a flxed value 
of time, the transit time resident pdf is constant. This shows that, in 1-D, 
the intensity of probability of finding water molecules at any position in the 
domain that transit at time t or less, given that t is fixed, is always identical. 
The transit time resident pdf gives the intensity of probability for the spatial 
position of water molecules, for a given value of transit time. In 1-D, this 
intensity of probability is uniform since the trajectory is unique. The resident 
age and life expectancy curves show the typical apparent discontinuity in 
concentration at inlet (resident pdf of A) and at outlet (resident pdf of E). 
These discontinuities have a drawback on the transit time resident pdf, which 
is not necessarily equal to the age resident pdf at outlet, and equivalently 
not necessarily equal to the life expectancy resident pdf at inlet. This can be 
attributed to dispersion effects at the boundaries. Since the Cauchy condition 
is homogeneous for T = 0+ at inlet for the age pdf problem, and at outlet 
for the life expectancy pdf problem, backward movement of water molecules 
by dispersion (i.e. in the reserved direction of velocity) is put down by non- 
zero age and life expectancy resident concentrations at the boundaries, the 
magnitude of which is higher the higher the dispersivity. If we consider e.g. 
the inlet boundary, the age and life expectancy concentrations are both not 
nil. The convolution of both pdfs is, therefore, not equal to any of them (age 
pdf at outlet, life expectancy pdf at inlet) since both concentrations can have 
a significant value at inlet and outlet at the same time. 

The spatial distribution of the transit time pdf is also ruled by a differen- 
tial equation. Combining Eqs. (2a) and (5a) after a Laplace transform, the 
following equation can also be obtained: 



where C^^ denotes the inverse Laplace transform. Eq. (14a) is the transit time 
pdf differential equation. It is of steady-state kind, with a source term that 
accounts for the divergence of the age and life expectancy dispersive fiuxes. For 




(14a) 



(14b) 
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the pure advective case in divergence- free flow fields, Eq. (14a) simplifies in 
c^-Vgrix, t) = 0. This local condition states that the flux vector and the transit 
time gradient have to be always perpendicular, as a requirement for keeping 
transit times constant along the flow paths. The resolution of Eq. (14a), which 
is of hyperbolic kind, is linked to technical difficulties, e.g. for the evaluation 
of the source term (14b) over the domain, and is not beneficial since the pdfs 
gA{x,t) and gE{x,t) must be evaluated in a preliminary step. 

2.5 Mean age, mean life expectancy and mean transit time. 

The average values of the probability density functions gA, gE and gr are 
defined by their first order temporal moment, the n^^ moment being 

M = u^guix, u)du = (-1)"^^^^^^^^ (15) 

with U = A, E, or T. Applying the convolution theorem in the Laplace domain 
results in the transformed Eq. (12), grix, s) — gA{x, s)gE{x, s). Accounting for 
the pdf property gu{x., s = 0) = 1 and enforcing Eq. (15) for n = 1 and 2 yields 
the average transit time and its variance: 

(T) = {A) + {E) (16) 

i^2[gT] = fi2[gA] + fi2[gE]+'2{A){E) (17) 

A9T] = c7^[9A] + (7'[gE] (18) 

with {A) = {A){x) = fiM, {E) = {E){x) = and (T) = {T){x) = 

liilgr]- The average age and average life expectancy of a water sample sum 
up to the average transit time. Since the variables A and E are supposed to 
be stochastically independent, the variances o"^ of the age and life expectancy 
pdfs also sum up to the variance of the transit time pdf. Using the operator 
in Eq. (15), Eqs. (2a) and (5a) can be transformed into their -n}^ normalized 
moment form. The recursive application of the operator (15) to Eqs. (2a) 
and (5 a) yields 

- V ■ CyXn[gA] + V • BVUnigA] - QofinigA] + ^n/^n-lifi'^] = (19) 

for the forward n*^ moment ADE, and 

V • (IfinlgE] + V • DV/inbij] - qi^n[gE] + ^^/^n-l bi?] = (20) 

for the backward n^^ moment ADE. The n^^ moment forms of the ADEs (2a) 
and (5a) are only dependent on the (n — 1)*^ moment For instance, since 
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A*o[fl'A] = 1 and gA{x,0) — 0, the first moment form of Eq. (2a) (for n = 1) 
corresponds to the mean age equation as defined by Goode [35], in which the 
mean age is the average over a sample of water molecules of the time elapsed 
since recharge: 

-V-q(A) + V-DV(A)-go(^)+0 = O in O (21) 

The first moment form of Eq. (5a) gives the backward adjoint mean life ex- 
pectancy equation: 

V -qiE) +V ■'DV{E) -qi{E) +(1)^0 in Q (22) 

Finally, the mean transit time equation is deduced by subtracting Eq. (21) 
and Eq. (22): 

ci-V{T) = {Sd) (23a) 
(S-rf) = V-DV(A) - V-DV(£;) (23b) 

where the divergence of the advection term has been developed in order to 
annihilate the fiuid source and sink terms. The boundary value problems (2) 
and (5) involve that the finite moments of the age and life expectancy pdfs exist 
for homogeneous boundary conditions. By definition, the mean groundwater 
age in a steady-flow reservoir, or mean residence time, can be determined from 
the average of the normalized flux concentration response to a narrow fiux in- 
put uniformly applied on the recharge limits, since this breakthrough curve 
corresponds to the water molecules residence time distribution [21,41,39]. The 
mean groundwater age can then be calculated by prescribing at all infiow 
boundaries a solute mass that is proportional to the water flux [37]. Conse- 
quently, Eqs. (21) and (22) can be solved by assigning (A) = on the inlet 
limits r_, and (E) = on the outlet limits r+, respectively. Mean age and 
mean life expectancy are continuously generated during groundwater flow, 
since porosity acts as a source term in Eqs. (21) and (22). This source term 
indicates that groundwater is aging one unit per unit time, in average. The 
mean age and mean life expectancy equations can be easily handled by numer- 
ical codes that solve ADEs, by distributing a source term equal to porosity, 
and by reversing the velocity field for the case of Eq. (22). Eq. (23a) would 
require the boundary conditions (T) = (E) on r_, and (T) = (A) on r+. 
However, mean transit time can rather be obtained by solving the Eqs. (21) 
and (22), and by post-processing Eq. (16). If dispersion is nil Eq. (23a) is 
simply q ■ V(T) = 0, the solution of which is comparable to the stream lines 
in a fiow model, and associates to each line the advective transit time from 
inlet to outlet. 

According to Parker and van Genuchten [53,54], Kreft and Zuber [41,42], and 
Sposito and Barry [65], temporal moments have a real physical meaning if 
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they are related to flux concentration fTinctions, while spatial moments must 
characterize resident concentration functions. Flux and resident concentration 
depend themselves on the measurement technique [65,59]. An age date as 
deduced from isotopic measurements corresponds to an age average over the 
number of molecules in the water sample, and may often be close to a resident 
concentration. Mean age computations using Eq. (21) are, therefore, well- 
suited for fitting isotopic age dates. 



2.6 Theoretical illustration example. 

The numerical solutions presented in this work have been performed using the 
Laplace Transform Galerkin (LTG) finite element technique [68], which allows 
eliminating the time-derivative in Eqs. (2a) and (5a). The Cauchy type bound- 
ary conditions (2c) and (5c) imply that Eqs. (2a) and (5a) must be formulated 
according to their divergence form [23] , such that a total mass fiux continuity 
at the boundaries must be properly handled. The assembled linear system of 
equations is solved by the accurate ILUT-based sparse iterative solver [62] 
with complex arguments and GMRES [61] or BiCGSTAB [70] acceleration. 
The numerical inversion of the Laplace transformed functions is performed by 
the algorithm of Crump [15]. The quotient- difference algorithm is used to ac- 
celerate the convergence of the Fourier series [24]. This algorithm proved to be 
very efficient to treat inversion at the neighborhood of discontinuities or sharp 
fronts, and the required computational effort, which is linearly proportional 
to the 2n + 1 number of discrete Laplace variables, is highly diminished with 
respect to other acceleration methods. 

We use here the simple example of a theoretical vertical aquifer section (see 
Fig. 3) to illustrate age, life expectancy and transit time computations. The 
configuration of the flow problem corresponds to the typical case known as 
'well-mixed reservoirs', that generate an exponential-like transit time distri- 
bution at outlet. This aspect will be discussed in more details in Section 4.2. 
The aquifer is homogeneous, and is uniformly recharged on top by a constant 
infiltration rate. The outlet, which could represent a trench, is simulated by 
means of a prescribed constant hydraulic head from the top to the bottom 
of the section. We solve the boundary value problems (2) and (5), as well as 
Eqs. (21) and (22). The model is discretized into lOO'OOO bilinear quadrangles 
of size Ax = 0.5 m and Az = 0.25 m. Fig. 3a shows the fiow field by means of a 
forward particle-tracking representation. Fig. 3b gives the backward particle- 
tracking solution, which represents the purely advective life expectancy spa- 
tial distribution. In Fig. 3c, we have plotted the solutions of Eqs. (21), (22) 
and (16). Mean age and mean life expectancy are in very good agreement with 
particle-tracking solutions. Isolines of mean age are horizontally distributed in 
relation to the flow configuration: velocity increases upstream to downstream 
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(increase of hydraulic gradient towards outlet), but also becomes more and 
more horizontally distributed. Travel paths to reach a same depth are longer 
downstream than upstream, but since velocity increases downstream, it cre- 
ates a compensating effect and the corresponding travel times finally become 
similar, in average. Mean life expectancy gradients are mainly horizontal, and 
gradually increasing upstream. Mean transit time provides a good represen- 
tation of the flow field; its distribution is very close to the fiow lines shown 
in Fig. 3a. However, the mean transit time solution owns the additional in- 
formation on the total time required by water molecules to travel from inlet 
to outlet. Finally, in Fig. 3d are given some observed pdfs (see Fig. 3a for 
the location of the observation points). Successive age pdfs along a horizontal 
profile confirm the vertical gradient of mean age shown in Fig. 3c, each age pdf 
at a same depth being very comparable to one other. Along a vertical profile, 
successive age pdfs are more and more spread in relation to the effect of longi- 
tudinal dispersion, which is bigger the longer the travel path. The behavior of 
the life expectancy pdfs is similar to that of age, but in the reversed direction 
of velocity. The transit time pdfs along a horizontal profile are very similar 
to each other; they are simply shifted along the axis of time. This reveals a 
same amount of dispersion affecting transport of water molecules along each 
travel path from inlet to outlet. Vertically, the functions grix, t) show more 
and more dispersion with increasing depth, in relation to the vertical mixing 
of water molecules. 



3 Generalized Reservoir Theory. 

Recent studies [28,29] proposed a direct approach to calculate outlet transit 
time distributions by applying the reservoir theory (RT) [27,7] to the average 
groundwater age field, resulting from a pure advective transport solution. In 
the following, we show how the RT can be generalized to advective-dispersive 
systems by manipulating the forward and backward transport equations. We 
first introduce some characteristic reservoir time probability distributions. 
Two equivalent mathematical formulations are then proposed, in order to 
evaluate the discharge zone transit time pdf as a function of the reservoir 
internal physical properties. 

3.1 Characteristic reservoir distributions. 

When a specific age distribution is assigned to each elementary water volume 
in the reservoir, the volume of mobile water can be classified in a cumulative 
manner with respect to the age occurrence in the reservoir. Let M^(t) be the 
cumulated amount of water molecules in the reservoir with an age inferior 
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Fig. 3. Age, life expectancy and transit time computations in a 2-D vertical theo- 
retical aquifer, (a) Geometry, boundary conditions and flow field representation by 
forward particle-tracking with marker-isochrones each 1 year; (b) Backward parti- 
cle-tracking with marker-isochrones each 1 year; (c) Spatial distribution of mean 
age (solid lines) , mean life expectancy (dashed lines) and mean transit time (dotted 
lines) with 1 year of time-increment; (d) Age, life expectancy and transit time pdfs 
at some observations points. Parameters: K = 10~^ m/s, (p = 0.35, ol = 0.1 m, 
ar = 0.001 m, = 0. 

or equal to a particular value t, and mA{t) be the corresponding probability 
function, or reservoir internal age cdf. The function mA{t) is the porous vol- 
ume normalized function (t) , which is evaluated by integrating over O the 
probability of finding water molecules with an age t or less, assuming that 
each molecule has entered the system on r_: 

rriAit) = = ^ / <l>GA{x,t)dn = / <p(fgA{x,T)dr)dn (24) 

Mo Mo Jn Mo Jn \J0 J 



with Mo being the total amount of mobile water (aquifer porous volume). The 
function mA{t) cumulates the probability of finding water molecules which 
have travelled until a position x before time t. The function M^(t) is zero at 
origin and tends towards the total porous volume Mq at infinity. The internal 
age frequency distribution function ipAit) is the pdf associated to m^(i:), and 
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from Eq. (10) it follows 

Mt) = '-^ = ^Jj,.Mdn (25) 

Note that Eq. (25) corresponds to the zeroth spatial moment of the age pdf, 
which is equivalent to the age mass [59]. The function ipA{t) gives the probabil- 
ity density of finding elements in Q that have reached the age t, and ipA{t)dt is 
the probability that A lies in the interval [t, t+dt]. Thus, the quantity MoV^a(^) 
represents the number of elements per unit time (or flow rate fraction) that 
are in the interval [t, t + dt], and is equivalent to the zeroth spatial moment of 
the age pdf gA{x, t). We similarly deflne the internal life expectancy cdf mE{t) 
and pdf ipEit)'- 

Mt)-'^-j^JjMi^ (27) 

The function niEit) cumulates the probability that a water molecule in the 
reservoir will reach the outlet before time t. 

We finally consider the internal distribution of the groundwater molecules 
transit time pdf Qt as deduced from the convolution integral in Eq. (12), to 
introduce the functions mxit) and ^(t) as the internal transit time cdf and 
pdf, respectively: 

= = /n =wjn* ill *(^- ^'*) f^**' 

The function Mxit) corresponds to the mass of mobile water in the reservoir 
having a transit time inferior or equal to t. The function ^(t) characterizes 
the probability density of finding water molecules within the reservoir that 
have a transit time inferior or equal to t, and the quantity MQ^{t)dt gives 
the amount of water molecules that travel through Q within the time interval 
[t,t + dt\. 

3.2 The transit time pdf of discharge and recharge zones. 

A reservoir discharge zone is a particular portion of finite size, which inter- 
cepts the groundwater molecules that contribute to the outflow rate. These 
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water molecules have contrasted arrival times that converging water fluxes 
mix together before flowing out. At a given position on a discharge boundary, 
the transit time pdf is the age pdf. To characterize the contribution of each 
age flux event on an outflow boundary in terms of transit time probability, 
it is convenient to average the age probability fluxes on r+. Under steady- 
flow regime, the representative transit time distribution 9?a(^) of the reservoir 
outlet zone can be defined as a fiux averaged concentration [59], i.e. v^a(^) is 
evaluated as the flow rate-normalized sum on r_|_ of the total age mass flux 
response function J a resulting from a unit flux impulse on r_: 

<fA{t)^^[ JA-ndr^^f [qgA-'DV9A]-ndr (30) 



where n is a normal unit vector pointing outward the boundary, and where 
Fq represents the total discharge flow rate throughout the bounded domain, 
which at steady-state is evaluated by: 




Note that for the sake of simplicity, internal sources and sinks are neglected, 
Qi — Qo — ^- While flux concentrations are deflned with respect to a control 
plane orthogonal to the velocity vector direction, the outlet transit time prob- 
ability function in Eq. (30) is defined by the projection of the total age flux on 
the arbitrary-shaped boundary r_|_. Inserting Eq. (7) into Eq. (30) produces 
the equivalent relation: 



1 



ndr 



(32) 



With Eq. (32), one can see that, if the velocity vector q points in the direction 
of the outward normal unit vector n, then the dispersive term inside the 
brackets vanishes, and the pdf (pA{t) is equal to the total steady flow rate- 
normalized sum of the flux- weighted age flux concentration pdfs on the outlet 
limit. This is always the case in one-dimension. The cross product term in 
Eq. (32) reveals also that since the outflow limit is of arbitrary shape, then 
when velocity does not point in the direction of n a dispersive correction term 
is required. This is related to the fact that flux concentrations are defined with 
respect to a control plane which is orthogonal to the velocity direction. 

The discharge fiow rate can also be described as a cumulative function of the 
transit time values. We introduce the function /^(t) as the transit time cdf 
of the reservoir outlet, i.e. /a(^) is the probability that the molecules fiow 
out with a transit time t or less, such that it corresponds to the normalized 
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cumulated outflow function FA{t): 

Fo Jo 



(33) 



Note that the function /a(^) can also be viewed as the integral on r+ of 
the total mass flux deduced from the resident concentration solutions of the 
ADE (2a) for a unit step- input of mass flux on r_. By analogy, the life ex- 
pectancy pdf and cdf of the reservoir inlet r_ can similarly be defined by 

ipE{t) = ^ f JE-ndT^-^f [qgE + 'DVgE]- ndT 

(34) 



1 

FqJv 



= l\E{T)dT (35) 

Fo Jo 

Eq. (13) indicates that the transit time flux pdf equals both the life ex- 
pectancy flux pdf at inlet, and the age flux pdf at outlet. As a matter of 
fact, assuming neither addition nor subtraction of mass during transport, 
steady-flow conditions and stationary state, then, theoretically, every prob- 
ability flux on r_|_ has an identical counterpart J^; on r_ and vice- versa. 
Therefore, Eqs. (30), (32) and (34) relate to the same and unique function, 
(pit) = (pA{t) = ^E{t), and thus f{t) = fA{t) = fE{t). In fact, each element 
added in the reservoir at the inlet must exit at some position Xo of the outlet 
sooner or later. Each element added at the position Xo at the outlet must 
travel the same distance upstream, and thus spend the same time-span within 
the reservoir, backward in time, before reaching a position somewhere at the 
inlet limit. 



• ndr 



At the reservoir outlet, the arrival times qa are distributed along the discharge 
boundary, implying mixing and superposition of the information carried by 
each breakthrough curve. Moreover, within the reservoir mixing processes can 
also be important, and the true minimum and maximum ages can be diluted. 
To characterize an outlet representative transit time distribution, we must 
ensure that the minimum and the maximum ages are captured. Technical 
problems can often occur when solving Eq. (30) or Eq. (34), because mass 
flux line/surface integration is required. If the outlet is of small size, then the 
capture of calculated breakthrough curves, or the identification of particle ar- 
rivals, reveals to be ill-posed, mainly because of the loss of information due to 
the mixing of converging fluxes. Hence, numerical methods will generally re- 
quire a high level of refinement in the neighborhood of these integration limits, 
which rapidly becomes a computational limiting factor. Because the transit 
time pdf (p{t) on the inlet limit is identical on the outlet limit, discretization 



21 



methods imply that the number of observation nodes should be the same at 
the inlet and at the outlet, in order to be able to recover the same break- 
through curves. In other words, the temporal resolution of the curve, when a 
counting of the individual arrival times is performed, is a direct function of 
the spatial refinement in the vicinity of exit zones. The same restriction affects 
other simulation methods, such as the random-walk procedure. 

In the following, we propose an alternative approach that is relaxed from 
the above-mentioned practical problems. Eqs. (2) and (5) are considered to 
simulate the age and life expectancy probability distributions in the reservoir 
fl. Integrating Eq. (2a) over Q, and making use of the divergence theorem 
(/f2 V • Fdn = /r F • ndr) results in 

[qgA-'DWgA]-ndT+^J^ct>gAdn = - [q^^ - DV^a] • ndF (36) 



Normalizing Eq. (36) by the steady flow rate Fq and accounting for Eq. (25) 
and (30), Eq. (36) can be turned into the following form: 

<PA{t) + rj^^ = -F /r f"^^^ ~ ■ ""^^ ^^^^ 



with the quantity tq being the turnover time commonly defined at steady-state 
as the ratio of porous volume to flow rate: 

To = — (38) 



Similarly, the integration of the ADE (5a) has the form 

d 



[ [qgE + DV^b] ■ndT+- f (PgEdQ = f [qgE - DV^b] • ndV (39) 



Normalizing Eq. (39) by Fq and accounting for Eq. (27) and (34) yields 

<PE{t)+ro^^^^ f [qgE-T>VgE]-ndr (40) 
ot Pq Jr_ 



The boundary integrals in the right-hand sides of Eqs. (37) and (40) can 
be simplified by accounting for the boundary conditions (2c) and (5c). For 
instance, inserting the boundary condition (2c) into Eq. (37) reduces the 
boundary integral to — Fo5(i), and inserting the boundary condition (5c) into 
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Eq. (40) reduces the boundary integral to Fo5{t), and Eqs. (37) and (40) be- 
come: 

Mt)+ro^^^S{t) (41a) 
Mt) + ro^^ ^ S{t) (41b) 

With Eq. (41a) wc have recovered the RT formulation of Eriksson [27], in which 
the outlet zone transit time pdf (pA{t) is proportional to the first derivative of 
the internal age pdf ■0a (^), thus characterizing the probability for the water 
molecules of being removed from Q per unit time. Eq. (41b) is an equivalent 
formulation which relates the recharge boundary life expectancy distribution 
to the internal life expectancy distribution. Since the functions v^a(^) and 
(pE{t) are equal, it follows from Eqs. (41a) and (41b) that ipAit) = i'Eit), 
which allows writing the following general RT formulation: 

m+ro^^5{t) (42) 



Eq. (42) generalizes the RT to advective-dispersive solute transport processes, 
and is vahd for both age and hfe expectancy. If dispersion is set to zero, 
the function ■0(i) can be evaluated by integration of the field gA{x,t) — 
S{t — {A){x)) as proposed by Etcheverry and Perrochet [28,29]. Similarly, 
the function ilj{t) can be evaluated by integration of the field gE{x,t) = 
5{t - {E){x)). Since ip{t) = ipA{t) = '^^(t), it follows from Eq. (9) that 
M{t) = M^(t) = ME{t). This points out the importance in allowing age 
and life expectancy dispersive fluxes crossing naturally the outlet and inlet 
boundary portions, because the use of the homogeneous Neumann condition 
on r+ and may lead to different results for ip{t) if contrasted boundary config- 
urations and flow conditions exist. The third-kind boundary conditions (2c) 
and (5c) are the most meaningful conditions for solving the age and life ex- 
pectancy problems. They ensure a pure total flux pulse input entering the 
system at inlet for the age problem, and at outlet for the life expectancy 
problem. These conditions become homogeneous for t > (zero flux), and 
do not allow backward mass losses by dispersion. This would not be the case 
when using a Dirichlet type condition, which may lead to incorrect solute mass 
balances. For the one-dimensional case, the use of the Dirichlet condition per- 
mits the simulation of the age and life expectancy pdfs, but directly for flux 
concentration pdfs (see Appendix A). 

The fundamental relation between the outlet transit time cdf f{t) and the 
internal age pdf '4>{t) given by the RT is obtained after integration of Eq. (42): 

Fo - F{t) = = (43) 
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or 



(44) 



Eq. (43) indicates that the outflow of water molecules that leave the system 
through r_|_ with an age older than t is balanced by the number of elements 
per unit time within Q that reach the age t, i.e. in the interval [t, t + dt\. The 
corresponding volume of groundwater reaching the age t in 17 is MQip{t)dt. 
In other words, the flow rate fraction of age t or less at the outlet is a func- 
tion of the spatial occurrence in the reservoir of water molecules of age t or 
less. Since the function f{t) is zero at the origin, it follows that the value of 
%l){t) at origin is the turnover rate ctq ('0(0) = (Tq — I/tq), independently of 
the level of dispersion. Since f{t) is a cumulative function, then ■il){t) must 
be monotonically decreasing and M(t) must be an increasing function with 
monotonically decreasing increments [27]. The pdf ^(t) is constant from zero 
to the minimum age tmin at outlet, with -0(0, .. . ,tmm) = co, which throws 
light on the fact that the probability per unit time of finding elements in Q, 
that have reached the age t < tmin is certain. From Eq. (43) it follows that 
M{t) has a constant derivative equal to Mo'(/'(0) = Fq until the minimum age 
tmin is reached at outlet. Note that the same considerations can be made for 
fit) = fEit), m = Mt) and M{t) = MEit). 

Fig. 4 illustrates in 1-D the outlet (or inlet) transit time pdf, the internal 
age (or life expectancy) pdf, and the internal transit time pdf resulting from 
the analytical resolution of Eqs. (25), (29) and (42) (see Appendix A). Under 
pure convective transport conditions (dashed lines in Fig. 4), the function ip{t) 
equals the piston- flow transit time pdf 5{t — tq), and the function i/j{t) is the 
Heaviside function H(ro — t)/To, such that ilj{t) = ip{0) = 1/tq until tmin = To, 
and ilj{t) = after tmin. The first temporal moment of the pdf Lp{t) (average 
transit time at outlet rt) is dispersion- independent, and equals the turnover 
time tq. The average internal transit time Tit and the average internal age 
Ti (first temporal moments of \l/(t) and ip(t)) are dispersion dependent. In- 
creasing longitudinal dispersion (low Pcclct numbers) generates short arrival 
times and tailing effects (see the variances of the pdfs in Appendix A), and 
thus old arrival times at the outlet as well as old ages within the domain, 
which are visible on the three functions (p{t), ■0(0 and ^(t) for a range of 
Peclet numbers. The function xjj{t) is constant from until the minimum age 
at outlet. Since the Cauchy type condition prevents backward losses by dis- 
persion at a: = 0, the value of ■0(t) at the origin is always 1/tq for any Pcclct 
number (Fig. 4b). From the spatial organization of age or life expectancy oc- 
currence it is possible to predict the transit time distribution of a reservoir 
outlet (or equivalently the life expectancy distribution of a reservoir inlet), 
without the need of 'counting' the arrivals of the water molecules at a bound- 
ary of finite size. Thus, the pdf Lp{t) defined in Eqs. (30) or (34) as a pure 
boundary property becomes a property of the reservoir internal structure and 
hydro-dispersive characteristics. The information that can be lost when ip{t) 
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Fig. 4. Reservoir theory pdfs for a 1-D semi-infinite flow domain, as a function of 
the Peclet number for Pe = 5,10,25,50,100,250 and 500. (a) Outlet transit time 
pdf; (b) Internal age (or life expectancy) pdf; (c) Internal transit time pdf. Time is 
normalized by the average turnover time tq = L/v, x hy a characteristic length L, 
and Pe = Lv/D. 

is directly evaluated at the reservoir exit zone (or inlet zone) is recovered with 
Eq. (42). A far more accurate evaluation of ip{t) is thus achieved, for which 
the main operation (domain integrals (25) and (27)) is not time-consuming 
and can easily be implemented for one-, two- and three-dimensional systems. 
A convenient way to compute Eq. (42) is to work in the Laplace space, since 
it allows handling all time-dependent quantities in a quasi-analytical way. 



3.3 Temporal moments of the reservoir theory probability density functions. 



A direct consequence of the RT is that the expected value of the mean residence 
time, i.e. the average transit time Tt at outlet or inlet, equals the reservoir mean 
turnover time tq. This property can be found by calculating the first temporal 
moment of the transit time pdf (p{t), and by making use of Eq. (44): 

f+oo f+oo f+oo A/[ 

n = / t^{t)dt = [1 - f{t)]dt = To / = ro = ^ (45) 

JO JO JO -To 

where use has been made of the pdf property in Eq. (10). Since the reservoir is 
considered under steady flow conditions, internal average time characteristics 
can be defined. The temporal moments of the global functions ^/'(t) and \E'(t) 
have physical significance since resident concentration has been integrated in 
space. The mean internal age and the mean internal life expectancy Tie are 
deduced from the first temporal moment of the function ip{t). Integrating by 
parts and making use of the relation (42) results in: 
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Since ipA{t) = = i^(t), we have set ria = Tic = ri. Using Eq. (16), the 

average internal transit time r^t is deduced from the first temporal moment of 
the function ^'(t), and is found to be equal to twice the value of t^: 



r+oo ][ r f+oo ^ r 

rit= / t^{t)dt =—/(/./ tgrix, t)dtdn = — / (l){T)dn 
Jo Mq Jn Jo Mq Jn 



Using Eq. (42), integrating by parts and inserting Eq. (46), the second tem- 
poral moment and the variance of the distribution (p{t) reveal themselves as 
being functions of Tq and Ti (or Tit) only: 

r+oo 

= / r(p{t)dt = 2ToTi = ToTit (48a) 
Jo 

a^[ip] = To(2Ti - To) = To(Tit - To) (48b) 



Finally note that by using the same technique than in Eq. (46) , one can show 
than the second moment of ijj{t) is a function of the third moment of ip{t), 

IJ,3[ip] = 3Toll2bP]- 



3.4 Accuracy of the RT approach. 



The accuracy of the RT method compared to a classical evaluation of the tran- 
sit time pdf ip{t) at the outlet limit is illustrated in Fig. 5, using a four-layered 
vertical aquifer. The outlet zone is of very small size (Fig. 5a), which forces the 
information on age to be highly mixed. As attested by Fig. 5b, the individual 
age pdfs at outlet can be of very different shape. The differences between the 
two evaluation methods are very important. The ffux-weighted sum of the age 
mass fluxes monitored at the outlet nodes suffers from a loss of information 
induced by the mixing of converging fluxes at the outlet surroundings. With 
the RT, this information is recovered, since we ensure that each contribution 
to the outflow rate is accounted for. According to Eq. (45), the average res- 
idence time Tt must equal the turnover time tq\ it is clear with this example 
that this property is not satisfled by straight application of Eq. (30), while 
the RT provides a very accurate solution (Fig. 5c). 
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Fig. 5. Theoretical four-layered aquifer illustrating the accuracy of the RT com- 
pared to the classical direct evaluation at the outlet of the function (p{t). (a) Model 
geometry, boundary conditions and head solution in meters; (b) Age pdfs moni- 
tored at outlet, from which a flux- weighted average evaluation of the transit time 
pdf is performed; (c) Outlet transit time pdf evaluated at outlet and using the RT. 
Transport parameters: ai = 2.5 m, ar = O-Ola^jDm = 0. 
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4 Direct evaluation of aquifer water volumes versus age, life ex- 
pectancy and transit time. 

The specific groundwater volumes related to a given range of ages or residence 
times are important quantities to consider when addressing aquifer manage- 
ment strategies. Assessing the long-term evolution of groundwater chemistry, 
or defining corrective measures aiming at restoring groundwater quality after 
a contamination event, indeed requires appropriate age and volume-related 
information. In this section, we complete the framework of the RT by analyz- 
ing the transit time cdf f{t) to show how complex aquifer porous volumes can 
directly be quantified as a function of age, life expectancy, and transit time. 

4-1 Characteristic groundwater volume functions. 

Provided the residence time distribution ip{t) and the total discharge Fq are 
known, the characterization of internal groundwater volumes can be done in a 
relatively straightforward manner by a simple analysis of the function F{t) as 
defined in Eq. (33), F{t) = Fq /q ip{u)du. In Fig. 6, the theoretical shape of the 
cumulated outflow function F{t) (or cumulated inflow function) is represented. 
The areas va, vt, and vq in Fig. 6a characterize different groundwater volumes 
that can be defined as functions of time t (Fig. 6b). The area VA{t) represents 
the total groundwater volume F in Q of age inferior or equal to t but that 
will experience a transit time superior to t: 



In Eq. (49), the pdf ilj{t) is expressed given the relationship (43) between ilj{t) 
and F{t). When the minimum transit time tmin is not nil, the function v^it) 
contains the volume of age inferior or equal to time tmm, which is the area 
VAi{t) — tuiiniFo — F{t)), and the volume of age inferior or equal to t and 
superior to time tmin, which is the area VA2{t) = {t — tmin)(-fo — F(t)). The 
area vt = vrit) is the volume of groundwater that flows through VL with a 
transit time t or less: 



with the area VQ{t) being the amount of exfiltrated water having travelled from 
inlet to outlet during an observation time period t, 



VAit) = Vn{A<ta.ndT >t} = t[FQ - F{t)] = Motip{t) 



(49) 



VT{t) = Vh {T < = tF{t) - vo{t) 



(50) 




(51) 
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and where use has been made of Eqs. (25) and (43) to express the function 
M{t). The quantity tF(t) in Eq. (50) is the total amount of groundwater water 
in fl that reaches the age t or less on r_|_, plus the volume flowing out with an 
age t or less. The amount of groundwater water vrit) is nil until the minimum 
transit time tmin, and reaches the total porous volume at the maximum transit 
time tmax- The function VT{t) contains the volume of age inferior or equal to 
time tmini which is the area VTi{t) = tmm-^(^) in Fig- 6a, and the volume of 
age superior to time tram, which is the area VT2{t) = {t — tmm)F(t) — VQ{t) in 
Fig. 6a. Note also that the total amount of groundwater of age tmin or less is 
given by the sum vai + f ti = ^min-^o- 

Since the amount of groundwater VT{t) is the internal volume that will leave 
the reservoir up to time t, it equals the internal transit time cumulative dis- 
tribution function, Vxit) = MT{t), and it represents a part of the function 
M{t). The complementary part is the amount of groundwater water of age t 
or less and of transit time superior to t, namely the function VA{t) defined in 
Eq. (49): 

M{t) = Vn{A<t}= VAit) + VT{t) (52) 

With Eq. (52), one can express VA{t) as the difference between two increasing 
functions that both tend to the porous volume Mq at infinity. The function 
VA{t) is thus zero at the origin and at infinity. Since VT{t) is zero until the 
minimum transit time tmin, 'i^A(^) must equal M{t) between and tmm- Dur- 
ing this time-span, which can be taken as the signature of badly recharged 
and/or advection-dominated systems for significant values of tmm, these two 
functions have a constant derivative (see Fig. 6b) equal to the steady fiow 
rate Fq = MoiplO). The behavior of this function (number of peaks, compared 
durations of increasing and decreasing parts) is instructive as to the volu- 
metric proportions of groundwater remaining a long time in the system, or 
flowing quickly to the outlet. The time after which the function VA{t) starts 
to decrease gives information on the importance of the water volumes in the 
aquifer with long or short transit times. If this time is relatively young, then 
the aquifer may present a good turn-over property, and vice-versa. 

Differentiating Eq. (52) with respect to time, and accounting for Eqs. (25), (29) 
and (42), yields 

^(t) = jm (53) 

This fundamental relation includes all the features of the RT in the most 
compact form. Compared to the standard rule (42), Eq. (53) is a great im- 
provement. It is simpler and may provide the transit time distribution ip{t) 
with much higher resolution and accuracy, in relation with the fact that no 
differentiation between (p{t) and ^(i) is required. 

Note that if F{t) characterizes the inlet cumulated inflow rate, i.e. F{t) ~ 
FofEit), then each aquifer porous volume deflned above as a function of age 
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Fig. 6. Theoretical cumulated outflow function, internal age pdf (scaled by the 
porous volume) and internal groundwater volume functions: (a) Cumulated outflow 
function F{t), with the indicated areas va, vt and vq representing characteristic 
internal groundwater volumes relative to a value of time; (b) Aquifer porous volumes 
M(t), VAit) and VT^t) as a function of age and transit time. 

becomes a function of life expectancy. By analyzing the system outlet transit 
time cdf, the quantification of complex groundwater volumes with respect to 
age, life expectancy, and transit time, is straightforward, and provides im- 
portant practical system insights. An obvious application is the groundwater 
resources protection and vulnerability assessments. The underground storage 
of toxic wastes also requires the knowledge of groundwater volumes that are 
potentially subject to contamination, and assumptions on the time spent by 
these volumes within the reservoir until they reach an outlet area. 



4-2 Characteristic times versus aquifer particularities. 



Following Bolin and Rhode [7] , some aquifer configuration particular cases can 
be drawn from the first moments of (f{t) and ipit)- 

The first considerations correspond to the case for which the average tran- 
sit time and the average internal age are identical, Tt = ti. Using Eqs. (45) 
and (46), the condition for having n = tq = t; is /q°° t[ip(t) — 4'(t)]dt = 0, for 
which a sufficient condition is (p{t) = il){t) = ;^exp(^— The exponential 
form of the transit time pdf is a direct consequence of the first-order linear 
differential equation (42): if one of the two pdfs ip{t) and tp{t) has an expo- 
nential form, then the other must be identical. This is the exponential model, 
often termed the well-mixed model, which is mathematically equivalent to the 
unit response function of a well-mixed reservoir. In chemical engineering, this 
model is used for reactors inside which the age distribution of the elements 
is uniform, i.e. there is a perfect mixing of the elements. Although the ex- 
ponential model is widely used by hydrogeologists to simulate isotopic data, 
its application in aquifer systems must be handled carefully, since it involves 
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a large number of poorly realistic assumptions on the aquifer structure and 
recharge conditions. Eriksson [26] interpreted the exponential distribution of 
ages in groundwater as a consequence of an exponential decrease of poros- 
ity and permeability with depth. Luther and Haitjema [44] argued that the 
conditions on the validity of the exponential residence time distribution in 
porous media (horizontal, un-stratified and homogeneous aquifer with respect 
to porosity 0, recharge rate /, and saturated thickness H), can be relaxed if 
the parameters 4>, I, and H vary in a piecewise constant way, such that the 
ratio (pH/ 1 remains constant throughout the domain. This ratio characterizes 
the system turnover time, which means that each water sample taken from 
the reservoir must lead to a mean age that equals the aquifer mean turnover 
time. In nature however, such system configurations and conditions on mean 
age may hardly be found. Etcheverry [30] showed that a simple linear vari- 
ation of the thickness H significantly infiuences the shape of the theoretical 
exponential residence time distribution. 

We now consider the case where Tt > n. This case corresponds to the situation 
for which only few water molecules leave the aquifer rapidly after having en- 
tered. Confined aquifer conditions and/or very distant recharge and discharge 
zones are typical settings leading to these features. For such a configuration, 
the outlet transit time pdf (p{t) is generally small or nil for young ages, attest- 
ing the existence of a minimum duration for travelling from inlet to outlet. 
The functions VA_{t) and M{t) tend to remain identical until the minimum 
transit time is reached. After this date, the function Vj^{t) should decrease 
rapidly, refiecting the fact that after the shortest travel distance from source 
to sink has been covered, the outfiow of older water molecules is concentrated 
over a short time-span, in relation to the relative uniformity of the travel dis- 
tances within the system. A narrow triangle-shaped function VAit) is typical 
of aquifers with significant minimum transit times and bad turnover property. 
Note also that from Eq. (46) one can see that for systems with large turnover 
time, or relatively low hydro-dispersive properties (i.e. when a[ip] is small), 
the average age in the reservoir may also be smaller than the turnover time, 
with the lowest possible value Ti — ^ foi" the piston-flow conflguration only 
{a[ip] = 0). 

Finally, the case < r; corresponds to the situation for which important 
amounts of water molecules enter the aquifer and flow out relatively rapidly, 
while sufficient amounts of water stay long enough to increase the value of Tj. 
For such a configuration, (p{0) must be bigger than ■0(0), and the two curves 
must both decrease and cross each other at a certain date, after which ■0(i) 
is higher than ^p{t). This situation may be encountered when the source and 
sink zones are close to each other, or when superficial recharge is uniformly 
distributed, such that the fraction of young water is important at outlet, but 
when the heterogeneity of the velocity distribution is such that long travel 
paths might lead to old ages within the domain. Karstic systems are typical 
media where the case Tt < Ti is encountered, in relation to the effect of the 
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high velocities in the karstic network, which can carry water molecules more 
or less independently of the surrounding low permeability matrix. Considering 
Eq. (46), the case Tt < can occur in systems with small turnover time, or 
relatively high hydro-dispersive properties (large cr[99]). 



5 Analysis of dispersion and aquifer geometry effects on age dis- 
tributions. 

In the following, 2-D theoretical examples are presented to illustrate the pro- 
posed methods. Analytical and numerical solutions are provided in order to 
test the sensitivity of the probability and density functions to the advective 
and dispersive parameters, as well as to the aquifer structure. 

5.1 Two-dimensional single flow system aquifer. 

A 2-D half-circular crown-shaped reservoir is used to simulate a single-flow 

system aquifer, homogeneous with respect to porosity and hydraulic con- 
ductivity K. A positive head difference is applied between the recharge and 
discharge areas (Fig. 7a). This geometry is well suited for the derivation of 
analytical solutions [29]. Due to homogeneity and symmetry, the flow lines 
remain parallel to each other, flow being one-dimensional along the flow line 
coordinate. The spatial distributions of mean age and mean life expectancy 
can thus be solved analytically (see Appendix B), and reveal to be symmetric 
(Fig. 7a), yielding a co-centrical distribution of mean transit time. The porous 
volumes va (= vai + '^Ai) and vt (= vti + vt2)-i quantified using the cumula- 
tive outflow function F{t), can easily be identified in Fig. 7 by using the mean 
age and the mean transit time isochrones. Note that the special case of this 
theoretical aquifer allows a good representation of the groundwater internal 
volumes with the mean values of the global distributions. However, the mean 
of the age, life expectancy, and transit time distributions would generally give 
a rough representation only, when e.g. significant dispersive processes take 
place. Etcheverry and Perrochet [29] analyzed the pure advective case and 
found that ip{t) is proportional to and that il){t) is proportional to ln(l/t) 
between the minimum and the maximum transit time (see Fig. 8a). The func- 
tion ^'(t) is a constant from the minimum to the maximum transit time, which 
indicates that the aquifer volumes ranging between a unit increase of transit 
time remain constant (the aquifer volumes between two iso-contours in Fig. 7b 
are all the same). The case = in Fig. 8a was analyzed by means of ana- 
lytical solutions (Appendix B). For this system, the pdfs v'(t), ip{t) and \I'(t) 
show generally moderate fluctuations with respect to longitudinal dispersion. 
The main effect of the ai coefficient can be seen in the increase of the spread- 
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(a) (b) 

Fig. 7. Analytical solutions for mean age, mean life expectancy and mean tran- 
sit time for the 2-D vertical half-circular aquifer (Appendix B): (a) Mean age 
(solid lines) and mean life expectancy (dashed lines) isochrones with 50 days in- 
crement; (b) Mean transit time isochrones with 50 days increment. Parameters: 
AH = H2 — Hi = 100 m, K = m/s, 4> = 0.2, inner radius ro = 250 m, outer 
radius R = 1000 m, minimum transit time tmin = 143 days, maximum transit time 
imax = 2285 days, turnover time tq = rt = 772.5 days. 



ing of the probability distributions. The outlet transit time pdf shows younger 
arrival times when ai increases, but also a longer tail. The functions ilj{t) and 
^(i) are affected the same way by longitudinal dispersion. The average internal 
age Ti and the average internal transit time r^t are dispersion dependent; they 
increase linearly with a^, indicating that dispersion leads to a higher average 
age of the system. Note that the second temporal moments of the internal 
age and transit time distributions increase proportionally with the square of 
ai, while the variance of (f{t) increases linearly with a^. Contrasted behav- 
ior appears when the coupled effect of longitudinal and lateral dispersivity is 
taken into account (Fig. 8b). The case 7^ was analyzed using numeri- 
cal solutions, by increasing the ratio a = a^/aT- The tailing effect decreases 
from a starting situation (a = 00 in Fig. 8b) with increasing values of lateral 
dispersivity (decrease of a ratio in Fig. 8b). Lateral dispersivity induces the 
mixing of ages, old water molecules moving laterally between the flow lines 
and replacing younger molecules, and vice-versa. The internal groundwater 
volume functions also attest of the above-mentioned dispersion effects. For 
example, the function VA{t) shows long tails for large values of a^, reinforc- 
ing the observed behavior at outlet of the transit time pdf for old ages. The 
internal increase of water volumes that require long transit times to exit the 
aquifer due to an increase of has of course its consequence at outlet with 
older arrival times. As with the pdfs v?(t), ip{t) and ^I/(t), this tailing effect 
decreases when lateral dispersion is added. 
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Time [days] Time [days] Time [days] 

Fig. 8. RT solutions for the 2-D vertical half-circular aquifer: (a) Analytical solutions 
(see Appendix B) for the pdfs (p{t), ip{t) and ^(t), and for the internal groundwater 
volume functions (in % of Mq) M(t), VA{t) and VT{t), as a function of the ai 
coefficient {ot = -Dm = 0); (b) Numerical solutions as a function of the ratio 
a = ai/aT, with ol = 50 m. 

5.2 Vertical multi-layer aquifer. 

In this second theoretical example, we consider a four-layered vertical aquifer 
system, as illustrated in Fig. 9a. From the top to the bottom of the model, 
the layers have decreasing thicknesses and increasing pore velocities. The do- 
main is discretized into 30 '000 homogeneous bilinear quadrangles. The total 
porous volume Mq is 65312.5 m^. A constant input flow rate enters the sys- 
tem along the left limit, using imposed fluxes of varying intensity proportional 
to the different layers hydraulic conductivities. The system turnover time is 
To = 162.267 days. The outlet is simulated at the top of the right side by a 
constant hydraulic head along a relative small zone of 15 m. The permeability- 
porosity couples have been set in a way that the pore velocity contrasts involve 
specific ages within each layer (velocities Vi in Fig. 9a), creating specific ages 
within each layer, as illustrated in Fig. 9b. The fact that the layers thick- 
ness diminishes with depth, while the influx intensity increases, is meant to 
create arrival time peaks at outlet of comparable orders of magnitude. The 
temporal moments and r^t are good indicators of the dynamics of the global 
system. Because they are volume-averaged quantities, their magnitude will 
depend on the water quantities of contrasted age and transit time, which are 
directly related to the flow and transport dynamics. For this example, T; and 
Tit show only small variations with respect to dispersion (Table 1), because 
even if velocities between layers are contrasted, flow in the system is generally 
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H= 150 m 




Fig. 9. 2-D vertical multi-layered aquifer: (a) Boundary conditions and Laplacian 
flow field in meters; (b) Mean age isochrones in days, illustrating the pore velocity 
contrast effect on age transport {ul = Ax = = 2.5 m, ax = 0, = 0). 
Average pore velocities (in m/day): vi = 4.0, V2 = 6.5, U3 = 10.5, ^4 = 23.5. 

rapid. Longitudinal dispersion has the effect of making the system older, by 
creating long tails that can be observed at the reservoir outlet, but also on 
the internal age and internal transit time distributions (Figs. 10a and 10b). 
Spreading in the direction of velocity is of course the main cause for this. 
Lateral dispersion presents the property of homogenizing the ages by mixing 
water molecules of different ages. For a given level of longitudinal dispersion, 
increasing the ratio a = ai/dT (Fig. 11) makes the system younger (tailing 
is lowered), with diminishing values of the mean internal age tj and the mean 
internal transit time th (Table 1). Longitudinal dispersion induces spreading 
of solute particles, and thus variability around the plume center of gravity. 
When lateral dispersion is added, the mixing of water molecules between flow 
lines homogenizes the ages. When longitudinal spreading is important, the 
exchange surface along the plume body is extended, and mixing can then be 
expected to be of high magnitude. The internal transit time pdf \l/(t) is of 
particular interest for the characterization of the internal organization of the 
flow dynamics. Its shape and particularly the number of modes it shows has a 
direct consequence on the shape of the outlet transit time pdf. If we take the 
example of the case a = 00 (a^ = 0) in Fig. 11, the function \l/(t) exhibits 
as many peaks as the outlet transit time pdf ip(t). However, the magnitude of 
these peaks are different for the two functions, and discern the transit time 
frequencies at outlet from those of the system pore volume. The peak of max- 
imum intensity is the first one for ip{t), and the third one for "^{t). If one 
looks at the time value corresponding to the first peak of <f{t), say ti, then 
inside the domain the density of probability that the water molecules have a 
transit time equal to ti is smaller than the same density of probability at the 
reservoir outlet. This is due to the fact that \l/(t) deals with volume-average 
probabilities while (p{t) deals with fiux-average probabilities. The four con- 
secutive peaks of the outlet transit time pdf correspond to the four families 
of water molecules, which transit within the four layers of the system, even 
though some parts of the curve attest to some mixing effects (the density of 
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Fig. 10. Calculated pdfs, cdfs and groundwater volumes for the 2-D vertical mul- 
ti-layered aquifer as a function of dispersion, for the ratio a = oj^/aT = 10: (a) 
Outlet transit time pdf; (b) Internal age pdf and internal transit time pdf; (c) 
Outlet cumulated outflow function; (d) Groundwater volumes versus time (in % of 
Mo). 

Table 1 

Dispersive parameters related to Figs. 10 and 11, and calculated pdf statistics (units 
in meters and days). 



Parameters 


Fig. 10 


Fig. 11 


a = ai/oiT 
ai 
ax 


10 10 10 
5.0 10.0 20.0 
0.5 1.0 2.0 


oo 50 10 
5.0 5.0 5.0 
0.0 0.1 0.5 


Statistics 






a[ip\ 


75.710 75.955 79.871 
98.902 98.915 101.010 


81.766 79.786 75.940 
102.156 100.905 98.916 



probability is not necessarily zero between two peaks). This simple theoretical 
example shows the complexity that is to be expected for the nature of the age 
and transit time distributions. It underlines how the significance of average 
age can lead to erroneous interpretations, regarding the system dynamics (e.g. 
inferences on average velocity) or hydrogeological problems related to risk and 
vulnerability assessment. Whenever the outlet transit time pdf is multi-modal, 
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Fig. 11. Calculated pdfs, cdfs and groundwater volumes for the 2-D vertical multi- 
-layered aquifer as a function of the ratio a = OL/ctT- (a) Outlet transit time pdf; 
(b) Internal age pdf and internal transit time pdf; (c) Outlet cumulated outflow 
function; (d) Groundwater volumes versus time (in % of Mq). 



then the average transit time can definitely not be used as a reliable quan- 
tity representative of the distribution since, as illustrated here, mean ages at 
outlet are mostly related to the lowest densities of probability (see Figs. 10 
and 11, Table 1). Once more this consideration points out to the interpretation 
of average age measurements, but also mean age simulations, which should be 
done very carefully. Generally, attention is given to hydrodynamic dispersion, 
heterogeneity, and long distance travelling induced mixing [35,71,72], which 
are factors that represent a source of uncertainty for the age dating meth- 
ods and environmental tracer data interpretations. Here one can see that the 
geological, the structural and the hydrauhc boundary configurations can by 
themselves be responsible of unrepresentative mean ages, even in advection- 
dominated transport regimes. Even if the entire groundwater age distribution 
cannot be measured in the field, mean age dates can be of great help for model 
calibration purposes. However, mean ages should a priori not represent an ab- 
solute simulation answer, and the knowledge of the entire age distribution 
should be of prime interest in most cases. 
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6 Summciry and conclusions. 



(1) The probability distributional evolution of groundwater age and life ex- 
pectancy has been simulated by forward and backward transient advection- 
dispersion type equations, according to proper boundary conditions. For a 
given position in space, the age and the life expectancy pdfs are complemen- 
tary distributions. Straight application of the principle of superposition to 
these distributions results in an integrated distribution of transit times, which 
characterizes the probability of having a given transit time from the recharge 
zone to the discharge zone, together with the quantities associated to that 
transit time. The transit time pdf is simply the convolution integral of the 
age pdf and the life expectancy pdf, and tells about the entire history of wa- 
ter molecules since recharge until discharge. Age, life expectancy, and transit 
time pdfs provide different kind of information, and each of them can reveal 
to be more advantageous than the other one depending on the hydrogcological 
insights to be provided. For instance, life expectancy and transit time distri- 
butions are well suited for vulnerability assessment problems (e.g. well-head 
protection, or underground storage of toxic waste), and they allow the map- 
ping of different regions within a recharge zone, in terms of residence time in 
the aquifer and associated properties. 

(2) The results of Eriksson [27] have been recovered by manipulating the 
advection-dispersion equation, extending thus the RT to systems with signif- 
icant dispersive components. In the classical RT neither advection nor dis- 
persion was considered. We have demonstrated here that the RT is still valid 
in systems with spatially varying velocity fields and non-negligible disper- 
sive/diffusive effects. The RT is a simple one-dimensional formulation, time 
being the only dependent variable. The outlet transit time distribution is de- 
rived from the internal age distribution, and therefore has a much more refined 
resolution. In so doing, mixing of converging flow patterns near the outlet is ig- 
nored, and the maximum transit time is never smaller than the maximum age 
in the reservoir. The RT formulation also applies to the internal distribution 
of life expectancy, and can be used to evaluate recharge zones life expectancy 
pdfs. It has also been shown that the RT can be expressed in a more com- 
pact form, which relates the outlet transit time distribution to its internal 
counterpart, the internal distribution of transit times from inlet to outlet. 

(3) From the reservoir characteristic distributions, fundamental additional 
transient information on water volumes and water fluxes can be gained. From 
the outlet transit time cdf, specific groundwater porous volumes can directly 
be identified and quantified with respect to age, life expectancy and tran- 
sit time. These functions can be very useful for aquifer characterization and 
intrinsic vulnerability assessments. They can be used to easily evaluate the 
magnitudes of young and old groundwater volumes in the aquifer. 
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(4) Using analytical and numerical solutions for theoretical aquifer configu- 
rations, some effects of macro-dispersion on simulated age and transit time 
distributions could be underlined. For instance, longitudinal dispersion can 
have a significant aging effect, while lateral dispersion rejuvenates the system 
through transverse mixing. Like in previous studies, it has been shown that 
the average age resulting from dating-methods, or direct simulations, can lead 
to erroneous interpretations. It is a well-known fact that dispersion can in- 
duce a high variability of the age distribution around the average. Here we 
have shown that not only the mixing processes could make the average age 
insignificant, but also the geological structure and the geometry of the fiow 
patterns. 

(5) The proposed methodology can equivalently be implemented in one-, two- 
and three-dimensions, and has considerable technical and numerical advan- 
tages, which may be pivotal in handling very large systems. In fact, when the 
outlet transit time pdf is defined by integrating all hydro-dispersive properties 
over the entire fiow field, the level of refinement required by a stable transport 
model is generally sufficient. The RT ensures that the minimum and the max- 
imum age in the reservoir are captured at the outlet, which is hardly the case 
with traditional methods, mainly because of the mixing of converging fluxes 
near the outlet. In the present work, the RT has been combined with ADEs, 
and the equations were solved using the LTG technique. However, no restric- 
tion appears for the use of the methodology if age and life expectancy pdfs 
are calculated by other transport models, such as random- walk simulators. 

(6) The presented models have been developed for the global aquifer system, 
regardless the number of individual inlet and outlet zones. At this stage, the 
RT has been rendered applicable to the whole system. In a subsequent article 
(this issue), we generalize the RT to any observation zone, and to systems 
with several inlets and outlets. 
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Appendix A: Reservoir Theory for a 1-D semi-infinite domain. 



The age resident pdf is found by solving the one-dimensional form of the age 
pdf ADE (2a) with the Cauchy type condition vgA{0,t) — DdgA{0,t)/dx — 
v5{t), and with a zero concentration gradient at a point at infinity. The cor- 
responding age flux pdf is the flux concentration g-^, which is deduced using 
the one-dimensional form of Eq. (7), g^ = 9a — These solutions are e.g. 

given in Jury and Roth [39]. Using the dimensionless variables 

X ^ V ^ vL 
X=- T^-t Pe= — 
L L D 

where L is a characteristic length defining the supposed outlet position, and 
where Pe is the Peclet number, the age resident and flux pdfs read 

T) = ^ exp _ ^ ^^p(xPe) erfc 

/ Pe(X-T)2\ ^ 

The life expectancy resident and flux pdfs can be deduced from (A.l) and (A.2) 
by substitution of X by 1 — X and X by 1, respectively: 




VPe ^ Pe(l - X - T) 

QeO^, T) = -j== exp 



VT V 4T ; 



(l-X)VPe / Pe(l-X-T)2\ ^^^^ 

The transit time flux pdf is calculated by straight application in the time 
domain of the convolution integral in Eq. (12): 



gU^,T)= r gi{X,T')gi{X,T-T')dT' = ^^exp^ ^^^^ 
Jo 2./7fT2 



2v/^T^ 



4T 



= ^{(1,T)=5^(0,T) (A.5) 

The transit time flux pdf in Eq. (A.5) could have been deduced by substitution 
of X by 1 in Eq. (A.2), or by substitution of X by in Eq. (A.4). The transit 
time resident pdf is found by convoluting (A.l) and (A.3) in the Laplace 
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domain: 



^t(X, s) = qaO^, sygEiX, s) = ^^-i-^ exp [1 - 7]) (A.6) 
where s is the Laplace variable, g is s-transformed state of the function g, and 



with 7 = y (1 + 4s/Pe). The inversion of Eq. (A.6) yields 
^t(X,T) = ^t(T) 

where use has been made of the shifting theorem 

£-M/(a« + &)} = ^/ 

(Jj \(Jj 

and of the following Laplace transform: 

I """^f'^J I = (1 + &c + 2bH) eMb't + be) erfc + bVi^ 

2bt f 

The internal age, internal life expectancy and internal transit time pdfs are 
obtained from Eqs. (25), (27) and (29): 



V^(T)= f\A{X,T)dX^ f\E{X,T)dX 
Jo Jo 



*(T) = /' griX, T)dX = ^t(T) (A.9) 
Jo 

In 1-D, the reservoir theory is equivalent to the mass conservation relation 
between flux and resident concentrations given by Jury and Roth [39]: 



dij(T) d /•! VPe / Pe(l - T) 



2" 



2v^T^ 



4T 



= 5^(1,T)=^^(0,T) (A.IO) 
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The dimensionless outlet transit time pdf is simply the age flux concentration 
at X = 1, or equivalently the life expectancy flux concentration at X = 0. The 
mean of the pdfs <^(T), V'(T), and \E'(T) read 

Tt - = 1 Ti - = 2 + ^ Tit - = = 2Ti (A.ll) 

and the spreading of these pdfs is measured by their variance: 

2r n 2 2r ,n (Pe+6)2 2r.n 2Pe+6 , 



Appendix B: Reservoir Theory in a single flow system. 



Consider the crown-shaped aquifer geometry in Fig. 7. Since the flow lines 
remain parallel to each other in the entire domain, we assume a 1-D transport 
process at each point a; = on the flow line coordinate. We follow the same 
procedure than in Appendix A, given that the porous volume is Mq = (f)T:{E? — 
rQ)/2, the steady flow rate Fq = KAHln(R/ro)/-K, the dispersion coefficient 
D{r) = aLv[r), and the pore velocity v{r) = q{r)/(f) = KAH/(f)nr. The age 
resident and flux pdfs are deduced from Eqs. (A.l) and (A. 2) by substituting x 
by 9r, and the life expectancy resident and flux pdfs are deduced by replacing 
X by (tt — 9)r. The transit time flux pdf is flnally obtained by replacing x by 
vrr in Eq. (A. 2). The average age, life expectancy and transit time are given 
by the flrst moment of the corresponding flux pdfs: 

(B.1) 

(T) (r) = (A) (9, r) + (E) (0, r) = ^ (B.3) 

Since on the outlet line {9 = tt) velocity points in the direction of the outward 
unit vector, the discharge zone transit time pdf (fi{t) can be evaluated by 
averaging the age flux pdfs on the outlet line, or equivalently by averaging the 
life expectancy flux pdfs on the inlet line {9 = 0). The internal age (or life 
expectancy) pdf ip{t), and the internal transit time pdf \l'(t) are calculated 
following the domain integration in Eqs. (25), (27) and (29), by integrating 
between the inner radius Tq and the outer radius R, and by integrating between 
and 71 (see Fig. 7). The average transit time Tt equals the turnover time Tq 
(see Eq. (45)): 
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The average internal age (or life expectancy) ti and internal transit time = 
2ri vary linearly with a^: 

= "'^"^^ ^ UKAHiR^ - rl) ~ "^^"^ + ^""^^^ ~ "o^^ ^^'^^ 

The variance of </?(t) can then be deduced from Eqs. (B.4) and (B.5) by en- 
forcing Eq. (46). 
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